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Motivated by the spate of recent experimental and theoretical interest in Mott insulating S = 1 
triangular lattice magnets, we consider a model S = 1 Hamiltonian on a triangular lattice interacting 
with rotationally symmetric biquadratic interactions. We show that the partition function of this 
model can be expressed in terms of configurations of three colors of tightly-packed, closed loops with 
non-negative weights, which allows for efficient quantum Monte Carlo sampling on large lattices. We 
find the ground state has spin nematic order, i.e. it spontaneously breaks spin rotation symmetry 
but preserves time reversal symmetry. We present accurate results for the parameters of the low 
energy field theory, as well as finite-temperature thermodynamic functions. 



I. INTRODUCTION 

Magnetism in Mott insulators has blossomed into an 
exciting frontier of quantum condensed matter physics 
on both the experimental and the theoretical front^. The 
magnetism in such materials is usually modeled by lat- 
tice spin Hamiltonians. These are the simplest many- 
body problems but yet they can be notoriously complex. 
In most cases, one must resort to uncontrolled approxi- 
mations to capture the rich landscape of emergent phe- 
nomena. Given the central role they play, it is of great 
importance to develop a collection of spin models that 
can be studied in the thermodynamic limit in an unbi- 
ased manner. These controlled results can serve as a 
benchmark for approximate methods and are sometimes 
necessary to study certain non-perturbative phenomena. 
While in one dimension, recent advances in numerical 
methods have allowed the study of almost any Hamilto- 
nian^, in two and higher dimension quantum Monte Carlo 
of a small set of "sign-problem" free models provide us 
with the only unbiased view of the quantum physics of 
spin models in the thermodynamic limit^. 

The motivation for the "sign-problem" free model we 
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FIG. 1. Loop configurations. The cartoons show how the 
color {x,y,z) of a site and its nearest neighbor evolve in a 
step of imaginary time and how this affects loop decomposi- 
tion with (a) off-diagonal operators, e.g. \xx){yy\ (b) diagonal 
operators, e.g. \xx){xx\, and (c) no operator. In the MC up- 
dates, updating the colors of the loops converts between (a) 
and (b). Dividing or joining loops of the same color converts 
between (b) and (c). In combination these updates are er- 
godic. The horizontal bars denote the action of one of the 
terms in the bond operator, Eq. (Isl. 



consider here comes from two recently studied materials 
with Ni^+ ions living on a triangular lattice, NiGa2S4^] 
and Ba3NiSb209^. Both materials are Mott insulators 
with the Ni^+ ions in a 5 = 1 state. In surprising dis- 
coveries, it was found that both materials have gapless 
ground states but do not realize the 120° magnetically 
order state predicted by a semi-classical analysis of the 
antiferromagnetic Heisenberg model. Despite further in- 
vestigations, the low temperature phase in these materi- 
als is still under debate^. Prompted by the experimental 
work, there has been a number of theoretical studies of 

= 1 models on the triangular lattice, using various ap- 
proximate methods and exact diagonalization on lattice 
with up to 21 sites. Depending on the model or method 
chosen, researchers have found magnetic states, proxim- 
ity to quan tum c riticalit>EI various spin nematic^^^^ and 
spin liquid^inil. One popular theoretical rationalization 
for the unusual experimental behavior is the presence of 
a strong biquadratic interactiorP^. 

Inspired by the extensive experimental and theoretical 
work we study a = 1 model on the triangular lattice 
with pure biquadratic exchange. Remarkably, we find 
that this Hamiltonian does not suffer from the sign prob- 
lem, and exploit this fact to study systems without any 
approximation on lattices with more than 10^ spins. Us- 
ing our numerical results we confirm that the biquadratic 
model has a spin nematic ground state, consistent with 
previous approximate studies and small cluster exa ct di- 
agonalization studies of the same 

HamiltoniarPGS. We 
then make accurate estimates for the low energy param- 
eters of this model, thanks to the large system sizes made 
accesible by our method. 



II. MODEL 

The model we are interested in can be expressed in 
terms of the spin-1 Pauli matrices S as, 

H = -KY,[S,-S,y (1) 

(ij) 

where the sum is taken over pairs of nearest neighbor 
sites of a triangular lattice. It is well known that the 
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pure biquadratic model has a staggered SU(3) symme- 
try on bipartitie lattices. On the non-bipartite triangu- 
lar lattice studied here it has only the symmetry of spin 
rotations. In the usual labeling of the eigenstates on a 
site, the states |1), |0) and |I) are eigenstates of Sz with 
eigenvalues 1, and —1. When written in the form above, 
the model, Eq. ([T]) appears to have off-diagonal matrix 
elements with both positive and negative signs, violat- 
ing the Marshall sign rule. The mappings to follow are 
most transparent in the orthonormal basis of eigenstates 
of Sx^ Sy and Sz with eigenvalue 0. These are. 
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|y) = -^(|l) + |l)) 
\z) = |0) 



(2) 



It is possible to verify by inspection that in terms of these 
basis states the biquadratic operator can be re-written 
simply as a projector. 
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where \S) = ^ i^i this basis it is transparent 

that the Hamiltonian, Eq. ([T]), satisfies the Marshall sign 
criteria, all off-diagonal matrix elements are negative or 
zero, even on non-bipartite lattices. We note in that al- 
though as mentioned above the model does not have all 
non-positive off diagonal matrix elements in the standard 
Sz basis, this situation can be remedied by the phase ro- 
tation |0) ^ z|0), in which case \S) = M±ll^idill^ then 
from Eq. ([3| it has all its matrix elements non-positive. 

The Marshall condition implies that our model Eq. ([T]) 
has a positive definite path integral, if we use the 
\x)^\y)^\z) basis. For concreteness we use the stochastic 
series expansion (SSE) method to change from Hamil- 
tonian to path integral at finite temperature T = 
in which Z = t^Tr(^^j^ (The mapping below 
goes through equally well in the usual Trotter path in- 
tegral approach). An SSE configuration is specified by 
assigning to each H in the trace one of the terms of 
Eq. (|3| that contributes to the Hamiltonian and its lo- 
cation on the lattice, e.g. \xx)ij{zz\ij (see Fig. [T]). The 
partition function is a sum over all SSE configurations, 
which at any given step in the evaluation of the trace 
is described by specifying the basis state the system is 
in at that time slice, i.e., by stating which one of the 
three colors, x, y or z, is on each lattice point. Colors 
can change through the action of the off-diagonal terms 
in Eq. ([3| in which two neighboring lattice points with 
the same color can both switch at the next time slice to 
a new color [see Fig. [ija)]. It can now be shown that 
every SSE configuration is identified with a unique loop 
decomposition: loops are constructed by starting at a 
point in space-time and moving along the world-line in 
the time direction until a bond operator is encountered 
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FIG. 2. Finite-temperature specific heat, Cv (circles) and 
uniform susceptibility, Xu , (squares) of the S = 1 biquadratic 
model, Eq. ([T]) on a 64 x 64 triangular lattice. The results are 
representative of the thermodynamic limit. The single spin 
Curie law S{S + 1)/3T is shown for comparison as a dashed 
line. We have plotted to make both data sets fit on the 

same scale. 



at which point one traverses to the neighboring site on 
the bond and switches the direction of motion until the 
loop closes; loops join space-time points with the same 
color and every space- time point is connected to only one 
loop (see Ref.l^ for related details). Putting the entire 
imaginary time history together, it is easy to see that 
every world-line configuration maps uniquely to a config- 
uration of closed tightly packed loops, each of which is 
assigned one of three colors, x^y or z. The weight of any 
loop configuration from the stochastic series expansion is 
then given by the simple formula. 
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where the sum over / is a sum over all closely packed loop 
configurations, ni is the number of operator insertions 
(shown in Fig. [l] as horizontal blue bars) and 3-^^ is an 
entropic factor that counts the number of ways colors can 
be assigned to the A/] loops in 



III. QUANTUM MONTE CARLO 

We implement the described mapping of Eq. ([T]) to 
a loop model and sampled the SSE configurations with 
Monte Carlo updates. Loop configurations can be up- 
dated by joining loops of the same color, dividing a loop 
into two smaller loops and changing the color of a loop 
(see the caption of Fig.[T]). Such updates are ergodic and 
can be chosen with the proper weights to satisfy detailed 
balance. We note here that past quantum Monte-Carlo 
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FIG. 3. Equal time spin structure factor, (left) Re[5's(^] 
shown in the full first Brillouin zone of the triangular lattice 
on a 32 X 32 system at /3 = 32. (right) the peak height of the 
structure factor at the corner of the BZ (so-called "K" -point) 
as a function of inverse system size. The solid circles are the 
QMC data and the solid line is a numerical fit that shows the 
spin order vanishes in the thermodynamic limit. 



work has looked at the biquadratic interaction using a 
different algorithm in which S = 1 operators are decom- 
posed into two S = 1/2 operators as well as a loop al- 
gorithm similiar to the one described here, but in the 
basis. These studies were however restricted to bi- 
partite lattices (1-d chains, 2-d square and 3-d cubic lat- 
tices j™. 

One can now derive expressions in the loop language 
for various quantities of interest in the spin model. The 
total energy E and the specific heat Cy = can 
be measured in the usual manner of the SSE^^. The 



uniform susceptibility Xu 



can be 



shown to be exactly | {W[) y where W[ is the tem- 
poral winding number and i is the index which sums over 
the loops identified with a given configuration in space- 
time. The spin stiffness is related to the well known wind- 
ing number estimator ps = T (^f where is 
the spatial winding number. 

We begin by studying the specific heat and the suscep- 
tibility of our magnet as it is cooled to its ground state. 
The data in Fig.[2]on a 64 x 64 system is representative of 
the thermodynamic limit. The peak in Cy{T) at T ^ 0.5 
signals the lifting of the extensive entropy of the spins in 
the high-temperature phase. The observed simultaneous 
saturation of Xu (T) is the classic behavior expected of an 
anti-ferromagnet. Even so, a study of the spin structure 
factor Ss{q) = J^r ^^^'^ {S{r) • S{0))^ shows no peaks 
that scale to a finite value in the thermodynamic limit, 
indicating the absence of simple magnetic long range or- 
der as shown in Fig. [3) 

A natural guess for a magnetic state with no Bragg 
peaks in the spin structure factor is a spin nematic. Such 
a state breaks spin rotational symmetry but preserves 
time reversal symmetry. The order parameter for the 
spin nematic is a traceless symmetric matrix. 
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FIG. 4. Finite-size scaling at fixed Kf3/L — 1. The upper 
panel shows finite size scaling of the data for the nematic order 
parameter, which in the thermodynamic limit extrapolates to 
Oq = 0.1376(2). The lower panel shows the extrapolation of 
the spin stiffness to a value ps = 1.298(2). 



where are components of the S = 1 operator. Since 
the order parameter is even in the spin operators, its 
condensation does not imply breaking of time reversal. 
In order to test for this order we measure the corre- 
lation function Cgir) = ^^{Qaa{r)Qaa{0)), which is 
easy to measure in our basis since it is diagonal. It 
is possible to show that the 0(3) invariant correlation 
function is simply related by a constant, i.e., CQ{r) = 
1 Ec./3(Q«/3(r)Q/3a(0)). Fig. Qa) shows a plot of = 
-p- CQ{r) versus 1/L, keeping the ratio of KP/ L = 1 
fixed so that we probe the ground state behavior. The fi- 
nite value for Oq in the T = and thermodynamic limit 
(combined with no condensation in the magnetic struc- 
ture factor) proves unambiguously that the ground state 
of our model, Eq. Q, is a spin nematic. It is interesting 
to ask how much our ground state deviates from a simple 
product nematic state with no quantum fiuctuations, a 
wavefunction for which is simply, |^prod) = Hi 1^^ = 0)*- 
In this state the (Qa/?) = ^ais/^ — ^az^^z^ from which we 
conclude that (0^'°"^)^ = 4/15 ^ 0.26666.... This is the 
maximum value the order parameter can take, if there 
are no fluctuations. From the extrapolation in Fig. [4] we 
find that (Og)^ = 0.1376(2) in the real ground state, 
the reduction of about 50% being due to the effect of 
quantum fluctuations. 

The breaking of the continuous 0(3) spin rotational 
symmetry results in gapless spin waves. Analogous to 
past work on the Neel state^ the low energy effective 
theory of the Goldstone bosons at leading order is char- 
acterized completely by just two parameters, the spin 
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FIG. 5. Determination of the spin wave velocity, c. In 
the main panel we show data for the (VKJ) and (VKr ) on an 
L = 32 system as a function of /3/L. In the inset we show 
how the /3/L value at the crossing of the types of plots in 
the main panel scales with system size. These are finite-size 
estimates for 1/c, from which we obtain c = 1.869(4) in the 
thermodynamic limit. 

stiffness, ps, and the spin wave velocity, c. The stiffness 
can be estimated by the usual winding number estima- 
tor. Fig. [4] (b) shows a plot of ps as a function of 
at fixed Kp/L = 1. A finite thermodynamic value for ps 
provides further evidence for gapless spinful excitations 
at long wavelengths. The spin wave velocity is a number 
that converts the units of space into the units of time. 
An estimator for the spin wave velocity, c, is hence the 
ratio of L/P at which the anisotropic 2+1 dimensional 
system behaves cubic. We compute this ratio for each 
L by adjusting p until (W^^) = (V^r^)- Then we carry 



out a thermodynamic extrapolation of the finite size 
data to obtain an estimate for c. The data for the wind- 
ing numbers for a 32 x 32 system is shown as an example 
in main panel of Fig. [5] We fit polynomials through this 
data and estimate the location of the crossing point for 
each L studied. The value of L/KjS at the crossing is 
then plotted versus l/L in the inset. Extrapolating to 
the thermodynamic limit we obtain that the spin wave 
velocity, c = 1.869(4). The estimates of the order param- 
eter Qa(3^ the spin stiffness ps and the spin wave velocity 
c, give a complete description of the low energy effective 
field theory to leading order. 



IV. SUMMARY 

In conclusion, we have provided a model S = 1 hamil- 
tonian with biquadratic interactions that is sign-problem 
free on non-bipartite lattices. We have shown how the 
partition function of this model maps to a remarkably 
simple loop model with positive weights, see Eq. Q. In 
this publications we studied how this mapping has al- 
lowed for a detailed study of the spin nematic ground 
state on the triangular lattice relevant to many recent 
experimental and theoretical works and provided high 
precision estimates for the ground state paramaters. Ex- 
tensions of the current work to carry out unbiased numer- 
ical studies of the role of disorder in quantum spin ne- 
matics and deconfined critical points out of the nematic 
phase^^ are exciting directions for future research. An- 
other field of experiments in which the nematic phase for 
5* = 1 spins has been discussed is in ultra-cold atoms in 
optical lattices; it is hoped that our high precision re- 
sults may be useful for guiding and testing experimental 
efforts in this context^^. 
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